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ABSTRACT 

We perform a combined fit to angular power spectra of unresolved infrared (IR) point sources from 
the Planck satellite (at 217, 353, 545 and 857 GHz, over angular scales 100 < I < 2200), the Balloon- 
borne Large-Aperture Submillimeter Telescope (BLAST; 250, 350 and 500 /an; 1000 < I < 9000), and 
from correlating BLAST and Atacama Cosmology Telescope (ACT; 148 and 218 GHz) maps. We find 
that the clustered power over the range of angular scales and frequencies considered is well fit by a 
simple power law of the form Cf ust oc l~ n with n = 1.25 ± 0.06. While the IR sources are understood 
to lie at a range of redshifts, with a variety of dust properties, we find that the frequency dependence 
of the clustering power can be described by the square of a modified blackbody, v@B(y, T c s), with 
a single emissivity index (3 = 2.20 ± 0.07 and effective temperature T cf f = 9.7 K. Our predictions 
for the clustering amplitude are consistent with existing ACT and South Pole Telescope results at 
around 150 and 220 GHz, as is our prediction for the effective dust spectral index, which we find to 
be «i5o-220 = 3.68 ± 0.07 between 150 and 220 GHz. Our constraints on the clustering shape and 
frequency dependence can be used to model the IR clustering as a contaminant in Cosmic Microwave 
Background anisotropy measurements. The combined Planck and BLAST data also rule out a linear 
bias clustering model. 

Subject headings: cosmic background radiation - cosmology: observations - infrared: diffuse back- 
ground - infrared: galaxies - submillimeter: diffuse background 



1. INTRODUCTION 

The angular power spectrum of Cosmic Microwave 
Background (CMB) temperature fluctuations currently 
provides vital constraints on cosmologic al models (e.g., 
iKomatsu et alJ 1201 It lLarson et al.1 120 111 ) . Experiments 
including the Atacama Cosmology Telescope (ACT), 
South Pole Telescope (SPT), and Planck satellite are 
now probing the CMB temperatur e power spectrum on 
arcminute scales dDas et al.1 120111 : iKeisler et al.1 120111 : 
IPlanck Collaboration et al.1 120111 ) An impr oved mea- 
surement of the Silk damping tail ( Silk 1968) improves 
constraints on, for instance, the scale dependence of pri- 
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mordial fluctuations, important for testing inflationary 
models, the number o f relativistic species, and early- 
universe exotica (e.g., IKomatsu et al.l 120 111) . On ar- 
cminute scales the contribution to the angular power 
spectrum from the primary CMB fluctuations becomes 
subdominant to extragalactic foreg rounds including in- 
frare d and radio point s ources (e.g,. IWhite fc Maiumdarl 
120041: IRighi et al.ll2008D . and th e thermal and kinetic 
Suny aev Zel'dovich effects (SZ; iSunvaev fc Zel'dovichl 
1970). Extracting the CMB signal requires understand- 
ing how the contribution from these components varies 
with frequency and angular scale. 

The infrared (IR) point source foreground is un- 
derstood to originate from high-redshift (z ~ 1 ~ 4) 
star-forming galaxies whose rest-frame emission peaks 
in the far-infrared due to thermal e mission from dust 
grains illuminated by starlight (e.g.. iBond et al.l IT9861 
Il99l IHuehes et al1[l998l : IBlain et aLlfl999t IDraindl2003ir . 
While our work concerns observations made in the mm 
and sub-mm we refer to these dusty sources through- 
out as 'IR' sources. Thermal dust emission from star- 
forming galaxies is also an important comp onent of the 
cosm ic infrared background (CIB - e.g., iPuget et al.1 
119961 ). Galaxies trace the large-scale structure and so 
are clustered, with a scale-dependent contribution to 
the power spe ctrum in addition to Poisson shot-noise 
(|Peeblesl fl98(il) . C lustering of IR sources was there- 
fore expected (e.g.. iBondl 119961 and references therein; 
iNegrello et al.l 120071 ) and has been detected in Sp itzer 
Space Telescope data at 160 /im (jLagache et al.ll2007D . by 
the Balloon-borne Large- Aperture Submillimeter Tele- 
scope (BLAST) and Herschel Space Observatory at 
250 350 and 500 dViero et all [20ul iCoorav et all 
[2010: lAmblard et all |2011|), in the microwave sky at 
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around 150 and 220 GHz by SPT and ACT Eall et al.l 
[20101 IDunklev et all [Mil IShirokoff et all 120 111), and 



in ea rly data from Planck (|Planck Collaboration et all 
120111 hereafter Pll). Correlations between cluster- 
ing at differ e nt fre quencies have also been detected: 
lHaiian et al.l (|2012l . hereafter H12) measure significant 
levels of correlation between BLAST maps and ACT 
maps at 148 and 218 GHz, detecting a clustered com- 
ponent at 4cr. Pll also find a significant correlation be- 
tween the Planck 217 GHz maps and those at higher fre- 
quencies (353, 545 and 857 GHz). Studying the clustered 
power of dusty galaxies in the sub-mm is simpler than 
at the mm CMB bands because they are the dominant 
extragalactic signal. Knowing that significant overlap 
exists between the sub-mm galaxy population and those 
sources responsible for the IR foreground in the CMB 
bands suggests that information about the former can 
help our understanding of the latter. 

In this work we combine large- and small-scale power 
spectra from Planck, BLAST, and correlations between 
BLAST and ACT to estimate the amplitude and scale 
dependence of the angular power spectrum of clustered 
IR sources. We present a simple power-law template 
to model the clustered source contribution that may be 
marginalized over when estimating cosmological parame- 
ters from foreground - conta m inated CMB m a ps as in, for 
instance, IHall et all ((20Trl . IDunklev et all poll , and 
IKeisler et al.l (|2011h . 

Previous ACT and SPT results (|Dunklev et al.lf20Tlt 
iShirokoff et al.ll2011h have found that the parameters ex- 
tracted from their CMB spectra are not particularly de- 
pendent o n the model adopted for the IR clustered power 
(see also iSehgal et all [20101: iFowler et all l20lo]) . The 
ACT and SPT data sets are not yet complete; the final 
data will include more sky coverage as well as measure- 
ments from additional frequency channels. iMillea et al.l 
(|2012f ) find that modeling the IR point source cluster- 
ing incorrectly for the final combined Planck and ACT / 
SPT data sets could introduce a significant bias in cos- 
mological parameters (they estimate la based on the dis- 
crepancy between two different IR clustering models). It 
is therefore important to understand the scale and fre- 
quency dependence of the clustered power in preparation 
for this future analysis. Improving constraints on IR clus- 
tering will also help constrain the SZ power spectrum. 

In this work we are primarily concerned with the IR 
sources as a contaminant in CMB maps. While the phys- 
ical properties of this high-redshift star-forming popula- 
tion are important for understanding the star formation 
history and galaxy evolution, in this analysis we do not 
attempt to extract information about, for example, the 
rcdshift distribution of the sources or the dark matter 
halos they occupy. 

In Section 2 we describe the data we use for our fitting, 
in Section 3 we explain our assumptions and methods, 
results are presented in Section 4 and a conclusion follows 
in Section 5. 

2. DATA 

Throughout this work we use 'auto-spectrum' to refer 
to a power spectrum calculated by correlating two maps 
at the same frequency, and 'cross-spectrum' to refer to a 
spectrum calculated by cross-correlating maps at differ- 
ent frequencies. 'BLAST x ACT' is to be understood to 



refer to the BLAST / ACT cross-spectra and so on. 

We use the BLAST 250, 350 and 500 auto-spectra, 
250 x 350, 250 x 500 and 350 x 500 (im cross-spectra, 
and BLAST 250, 350 and 500 x ACT 148 and 
218 GHz cross-spectra from H12, and Planck 857, 545, 
353 and 217 GHz auto-spectra from Pll to construct 
the template. The ACT data are from the 2008 observ- 
ing season and the BLAST x ACT spectra were calcu- 
lated from the ~8.6 deg 2 common to both sets of maps, 
as described in HI 2. We take the quadrature sum of 
the statistical and beam systematic uncertainties in Ta- 
ble 4 of Pll as the error on each Planck data point, 
neglecting any possible correlation in the beam uncer- 
tainty across different angular scales. We find that al- 
lowing such a correlation has minimal effect on our re- 
sults (Section 4.1.3). The BLAST and ACT beam un- 
certainties are subdominant to the statistical (noise) un- 
certainty across the range of angular scales covered by 
the BLAST and BLAST x ACT data and so we like- 
wise neglect any correlation in the errors on these spec- 
tra. The temperature to flux conversion factors given 
in Table 4 of Pll assume a source SED that varies as 
I(v) oc v \ a correction is then applied to convert to 
the real flux units (see Sections 5. 5 of Pll and 7.4.2 of 
iPlanck HFI Core Team d~al1l2rjll"l) . 

We subtract the Galactic dust emission (cirrus) com- 
ponent in the BLAST and BLAST x ACT spectra as de- 
scribed in Section 4.3 of H12, assuming the cirrus contri- 
bution varies with angular scale as £~ 2 . These spectra 
are not very sensitive to the details of the cirrus treat- 
ment, since the y are from a rel atively cirrus- free patch 
of sky (see also iDas et all 120111) . We assume that any 
contribution to the power spectra other than that from 
IR galaxies (for instance, radio galaxy-IR or SZ-IR cor- 
relations) is negligible. Removal of power in the Planck 
spectra that is not from extragalactic IR point sources is 
described in Sections 2 and 3 of Pll. 

The Planck and BLAST auto-spectra are presented 
in Figure 1. We also show th e 218 GHz ACT power 
spectrum from IDas et al.1 (120111) and t he 220 GHz SPT 
spectrum from IShirokoff et al.l (|2011l ). A WMAP-7 
best- fit ACDM CMB power spectrum (jKomatsu et al.1 
120111) has been subtracted from these points; we have 
not subtracted any power from radio point sources or 
the kinetic Sunyaev-Zel'dovich effe ct, as these contribu- 
tions are likely to be subdominant (|Dunklev et al.ll2011l : 
IShirokoff et al.l [20 111 ). All the spectra show similar an- 
gular scale dependence despite spanning a broad range 
of frequencies and almost two decades in angular scale. 
Note that a color correction to account for the differ- 
ent bandpass filters is required to make the Planck and 
BLAST spectra directly comparable - see Section 4 for 
details. 

The Planck, BLAST and ACT maps are calibrated 
using comparisons to various measurements including 
the orbital CMB dipole, FIRAS data and planetary 
temperature. Uncertainty in these calibrations must 
be accounted for in order to perform joint fits to 
the power spectra extracted from the different maps. 
The absolute photometric calibration uncertainties are 
7% for the Planck 857 and 545 GHz maps, 2% for 
the Planck 353 and 217 GHz m aps (see Pll and 
IPlanck HFI Core Team et~aT1 l20ll . 9.5, 8 7 and 9.2% 
for the BLAST 250, 350 and 500 /^m maps (|Truch et al.l 
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Fig. 1. — Planck (857, 545, 353 and 217 GHz) and BLAST (250, 350 and 500 /an) IR point source power spectra. The spectra include 
both shot-noise and clustered components. While these data span a broad range of frequency and angular scale there is a notable similarity 
in the angular scale dependenc e. The Planck error bars include both the statistical uncertainties and estimates of the systematic b eam 
uncertainty given in Table 4 oflPlanck Collaboration et al.l 120111) . Data points at I > 2000 from ACT at 218 GHz IDas et al.1120111) and 
SPT at 220 GHz UShirokoff et a.1.1120111) have been included for comparison with the Planck 217 GHz spectrum. We have subtracted a 
WMAP-7 best-fit ACDM CMB component from these spectra. No corrections for the different bandpass filter profiles have been applied; 
when the different filters and photometric calibration uncertainties are accounted for the Planck and BLAST data are in good agreement 
(Section 4). 



l200l . 7% for t he ACT 218 GHz m ap and 2% for the ACT 
148 GHz map (|Haiian et al.ll20fo() . Furthermore, the un- 
certainties in the BLAST calibrations are highly corre- 
lated; this is because the dominant source of uncertainty 
is the spectral energy distribution (S EP) of the star use d 
for the calibration of all three bands (jTruch et alll2009l ). 

3. POWER-LAW CLUSTERING TEMPLATE 

In t he halo model formalism (e.g.. IPeacock fc Smith! 
l2000l IScoccimarro et al.l [200lt ICoorav fc Shethl 12003 
see also earlier work, e.g.. IBondl 119961 and references 
therein), the two-point function of galaxy clustering is 
dominated on large scales by pairs of sources in different 
dark matter halos (the 'two-halo' term - roughly corre- 
sponding to the linear clustering regime), while on small 
scales galaxy pairs occupying the same halo ('one-halo' 
term - non- linear regime) become dominant. Although 
these two components have different scale dependence, 
an angular cor relation functi on varying as uu(9) oc 8~ d 
with (5 ~ 0.8 ((Peebles! 119801 corresponding to cluster- 
ing power Cf ust oc £~ 1 - 2 ) has been found to adequately 
describe the clus tering of, for instance , Lyman Break 
galaxies at z ~ 3 (iGiayalisco et alj|1998l) and local SDSS 
galaxies (jZehavi et al.ll2002f) . Significant deviations from 
power-law behavior have however been observed in anal- 
yses of more recent and comprehensive SDSS data (e.g., 



iZehavi et a fll200ll [20051 iBlake et alj|2008ft. as well as in 
high-redshift ga l axies (lOuchi et alJl2005HLee et alj|2006t 
iCoil et al.ll200a iWake et all 120 111) . For a recent discus- 
sion of the physic al origins of p o wer-la w galaxy correla- 
tion functions see I Watson et all (t201lh . 

Pll and lAmblard et al.l ()201lT ) find that a power law 
provides an adequate fit to the existing Planck and Her- 
schel /SPIRE unresolved IR point source spectra (see Ta- 
bles 6 and SI in those papers, respectively). This would 
suggest that IR spectra are not yet of sufficient quality 
to reveal deviations from power-law clustering behavior 
and we therefore also adopt a power law for the clustering 
component in this work. 

We model the total IR point source power spectrum 
from correlating maps at frequency v\ and V2 (zt = ^2 
for the auto-spectra) as 



(1/1,1/2) = A c (vi,v 2 ) 



+ CpK,z/ 2 ), (1) 



where £ is the multipole moment, A c and n are the clus- 
tering amplitude and index, Cp is the Poisson shot-noise 
and £q = 3000 i s the pivot scale. Th is differs from the 
form adopted bv lAmblard et al.l (|2011l) only in the choice 
of pivot £q. Unlike Pll we model the clustering and shot- 
noise as separate, independent components rather than 
modeling their sum as one power law. Motivated by the 
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apparent uniformity in angular scale dependence across 
the different spectra (see Figure 1) we fit for a single, 
frequency-independent, value of n. We fit for a separate 
shot-noise level in each of the 16 Planck, BLAST and 
BLAST x ACT spectra (see Section 3.2). 

The SED of the CIB, over the range of frequencies 
considered, has been found to be well-described by a 
modified blackbody of the form fe . g.. [Fixsen et al.lll998t 
lLagache et al.lll999tlGispert et al.ll2"000D : 

IcibO) otv B(v,T eS ). (2) 

In reality, the sources making up the CIB lie at a range 
of redshifts, wit h a variety of luminosities and dust tem- 
peratures (e.g . . | Haiman k, Knoxll2000t IKnox et al.ll2001l: 
Coppin et all 120081: iPascale et all 120091: iHwang et al.l 
2010D . Equation (2) is thus an approximation to the 
true CIB SED, which consists of the sum of many differ- 
ent (approximately) modified blackbody spectra, and the 
quantities /? and T e g are not to be interpreted as physical 
parameters. 

Pll found that the SED of the CIB anisotropics 
measured by Plan ck is also well-described by the 
iGispert et al.l (|2000l ) modified blackbody with emissiv- 
ity spectral index (3 = 1.4 ± 0.2 and effective tempera- 
ture T e g = 13.6 ±1.5 K. We assume the clustering power 
SED can also be described in this form, and parameterise 
the frequency dependence of the auto-spectrum cluster- 
ing power amplitude as A c = (I c ) 2 where 

/cM=/oQ B{u,T eS ), (3) 

with a single emissivity index and effective temperature. 

Different frequencies are sensitive to IR sources at 
different redshifts, with the importance of higher- 
redshift sources increasing at lower fre q uencies (e.g., 
Haiman fc Knoxl 120001 IKnox et al.1 120011 : iChapin et all 
2009HMarsden et al.H2009D . As a result we would not ex- 
pect there to be 100% correlation between the IR sources 
in different bands (particularly between the widely 
spaced ACT and BLAST bands) and we therefore include 
a clustering correlation factor f co „ as a free parameter 
for each cross-spectrum. The cross-spectrum clustering 
amplitude is then A c (vi,v 2 ) = /corral, v 2 )h{vi)I c {v2)- 

3.1. Shot-noise 

The auto-spectru m shot- noise is given by (e.g., 
IScott fc Whrt3[T999T) 

//■Scut d 2 /V 
dz j dSS 2 ^-(S,z), (4) 

where S cut is the flux cut applied to the map and 
d 2 N/dSdz are the differential source counts. 

In Pll the shot-noise levels were fixed using the IR 
galaxy evolution model of Bethcrmin et al. (2011 - here- 
after Bll). This model parameterises the evolution of the 
galaxy luminosity function and is fit to number counts 
over a wide range of IR wavelengths. While the model 
broadly fits the available data, there are discrepancies. 
The model undcrpredicts by ~40 % the shot-noise l evels 
measured at 500 u rn by BLAST (IViero et al.l 12009ft and 
220 GHz by SPT (jHa'll et al.l I2010D . We do not apply 
any priors on the shot-noise levels when constructing our 



template. We may, however, expect a considerable de- 
generacy between the shape of the clustering component 
and size of the shot-noise for the Planck data, because 
Planck is not able to probe the small scales where the 
shot-noise becomes dominant. As a result we consider 
the effect of adopting the Bll model predictions as pri- 
ors on the Planck shot-noise levels in Section 4.1.1. 

In terms of the source flux a nd S cu t we can write 
the clustering p ower as (e.g., iTegmark et all 12001 
IViero et al.ll2009l) 

Cf^ = /d*(^(*)) 1 (^ z {z)j P^k,z)\ k=l/X(z) , 

(5) 

where x is the comoving distance, dV/dz = x 2 dx/dz the 
comoving volume element, P ga i the IR galaxy power spec- 
trum, and the redshift-distribution of the flux, dS/dz, is 
given by 

— (z)= dSS——(S,z). (6) 

dz J Q dbdz 

The factor of S 2 in equation (4), compared to (dS/dz) 2 
in equation (5), suggests that the removal of the highest- 
flux sources will have much less impact on the clustered 
power than on the shot-noise. The Bll model predicts 
that, for instance, applying a flux cut of 250 mJy (the 
cut applied to the BLAST 350 /im map in H12) to the 
Planck 857 GHz map would result in a ^10% reduction 
in the shot-noise level compared to the Planck flux cut of 
710 mJy, but that (dS/dz) 2 would be reduced by < 1% 
at z ~ 0.2 and virtually unaffected for z > 1. We allow 
for the dependence of the shot-noise levels on flux cut 
by fitting separate shot-noise levels for each spectrum, 
as described in the next section. We assume that the ef- 
fect on the clustering power from applying different flux 
cuts is negligible for the data considered. Studies using 
higher-resolution Herschel maps will be able to investi- 
gate the dependence of the clustering power on flux cut 
in more detail. 

3.2. MCMC fitting 

We perform a simultaneous fit to the seven auto- 
spectra (four Planck and three BLAST) plus the nine 
cross-spectra (three BLAST x BLAST and six BLAST 
x ACT). The model spectra are binned for comparison 
to each data spectrum, with likelihood 

16 

-2lnL(d\6) = £[Cf ta - Cr dc \0)} 2 /cr 2 , (7) 

»=i 

for model parameters 9, data vector C data for the ith 
spectrum, and binned model spectra Cf lodel (9). Covari- 
ances between the sixteen spectra are neglected. 

We find that the clustering SED parameters I , (3 and 
T e g (equation 3) are highly degenerate. Only two of these 
parameters are really independent and so we fix T g to 
the best- fit value, 9.7 K. We choose vq = 530 GHz to 
minimise the degeneracy between Iq and /3. 

We therefore fit for 37 parameters: n, the clustering in- 
dex, Iq, /3, 16 shot-noise levels (one for each spectrum), 
nine cross-spectrum correlation factors (one for each 
cross-spectrum) and nine photometric calibration param- 
eters (one for each Planck, BLAST and ACT band). For 
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each calibration factor we enforce a Gaussian prior cen- 
tered at unity, with spread given by the nominal uncer- 
tainty listed in Section 2. The c ovariance be t ween the 
BLAST calibration factors from iTruch et all (|2009f ) is 
included. Apart from the calibration factors, all priors 
are uniform. 

Given the high dimensionality, we estimate the poste- 
rior probability distributi on using Markov Chain Monte 
Carlo (MCMC) analysis ^Metropolis et al.lll953l) . using 
th e sampling me t hods and convergence test described 
in iDunklev et al.1 (|2005[ ) . Chains used for analysis are 
about 10 6 steps in length, and are used to calculate one- 
dimensional marginalized parameter values and errors. 

To judge the goodness-of-fit of the model, we are also 
interested in the maximum likelihood. However, the peak 
of the likelihood distribution occupies only a small part 
of this high-dimensional parameter space, so the mini- 
mum \ 2 sampled in a chain of about a million steps is 
significantly larger than the true value (Ax 2 ~ 10 for a 
simulated 37-dimensional Gaussian distribution). There 
are numerous statistical methods to find the true peak 
of a distribution. We adopt a simple modification to the 
Metropolis algorithm in which chains start from the best- 
fitting point and then make a step in parameter space 
only when the posterior is improved, using a reduced trial 
step size. We have tested this prescription on simulated 
Gaussian distributions, and the peak is found to within 
A\ 2 — 0.1 in ~ 2 x 10 4 steps for 37 dimensions. We em- 
phasize that this modified code is used only to assess the 
global goodncss-of-fit, not to estimate the marginalized 
parameter distributions. 

To estimate the model spectra at each frequency, it 
is not sufficient to use the nominal values of v for the 
BLAST and ACT bands in equation (3), due to the fi- 
nite width of each filter. The correct effective v values 
depend on the SED of the emission mechanism. The 
clustering SED is estimated iterativcly by repeating the 
MCMC fitting; each time the best-fit clustering SED is 
integrated through the BLAST and ACT filter profiles 
and an improved estimate of the effective v values ob- 
tained. We found that after three iterations the SED 
converged, with effective frequencies of 1248, 829 and 
607 GHz (effective wavelengths 240, 362 and 494 fim) 
for the nominal BLAST 250, 350 and 500 bands, and 
220 and 150 GHz for the ACT 218 and 148 GHz bands. 
While the clustering SED is rising steeply at the ACT 
bands, the shift from nominal frequency is reduced by the 
narrowness of the ACT filters. For the Planck data the 
nominal frequencies are used, because the temperature- 
to-flux conversion factors provided in Pll are given such 
that the flux is correct at the nominal frequencies. 

4. RESULTS 

The model fits the data well. We find a best-fit x 2 °f 
132 for 122 degrees of freedom (150 data points minus 
28 parameters; we have not counted the calibration pa- 
rameters since they are strongly constrained by priors), 
giving reduccd-x 2 = 1.08. The marginalized mean val- 
ues of the clustered IR template parameters are given in 
Table 1. No errors are given for T g, vq or £q since these 
parameters are fixed as described in Section 3.2. Figure 
2 shows the best-fit clustered IR source power for each 
of the bands included in the fit. The best-fit shot-noise 
levels have been subtracted. The \ 2 contribution from 



each individual spectrum is also shown; in addition to 
the total x 2 being good we find that there are no spectra 
which are not individually well-fit by the model. 

We plot the BLAST 350 and 500 /xm spectra on the 
same axes as the Planck 857 and 545 GHz spectra. The 
BLAST spectra are color-corrected to account for the 
difference in the bandpass filters by taking the ratio of 
our clustering power predictions for the relevant Planck 
and BLAST band. We find these color correction factors 
to be 1.07 and 0.63 for the 350 and 500 (im BLAST 
spectra, respectively. Pll found values of 1. 05 and 0.7 by 
integrating the SED of IGispert et al.l (|2000[ ) through the 
BLAST and Planck filters. These values differ slightly 
because our best-fit clustering SED is slightly different - 
see Figure 3 for a comparison. 



TABLE 1 
Power-law clustering template 



Parameter 



Value 



1.25 ±0.06 
la (Jy sr" 1 ) 21 (2.43 ± 0.16) x 10~ 9 
2.20 ±0.07 
9.7 
530 
3000 



/3 

Teff (K) 

ko (GHz) 

e 



x'/d.o.f. 



132/122 



a A strong (89%) anti-correlation exists between n and Iq 

The clustering power at £ = 3000 is shown as a func- 
tion of freque n cy in Figure 3. Also shown is the scaled 
IGispert et all (|2000l ) SED (dashed line), which was fit 
to FIRAS and DIRBE measurements of the CIB inten- 
sity spectrum. It falls off more slowly with decreas- 
ing frequency than our clustering SED. Some differ- 
ence between the CIB mean and anisotropy SED may 
be expected, since the contribution of a source to the 
anisotropy SED depends on its clustering as well as spec- 
tral properties. 

IR clustering power predictions at £ = 2000 and 3000, 
calculated from our template for several Planck, ACT 
and SPT bands, are given in Table 2. The conversion 
from flux to temperature units was calculated in each 
case by integrating I c from equation (3) through the rele- 
vant bandpass filter. The effective frequencies (i.e., single 
frequency value that gives the same clustering amplitude 
as integrating through the filter) arc also given. Since 
the IR clustering amplitude is rising strongly with fre- 
quency at the CMB bands (see Figure 3), the different 
filter profiles lead to significant differences in the clus- 
tering power even for bands with closely-spaced nominal 
frequencies. The amplitude of the clustered power may 
vary by a factor of 10 across a single filter, as shown in 
Figure 3. 

The predictions at £ = 3000 can be compared to the 
ACT results from IDunklev et al.l (f2lml ) and SPT results 
from iShirokoff et al.l (|2011l) . ACT find a best-fit clus- 
tering amplitude of 4.6 ± 1.1 and 54 ± 13 piK 2 at 148 
and 218 GHz, respectively (we have calculated these un- 
certainties by adding their statistical and systematic er- 
ror estimates in quadrature). SPT find 6.1 ± 0.8 and 
57 ± 8 fiK 2 at 150 and 220 GHz, respectively, for their 
base-line model, which includes a power-law IR cluster- 
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Auto-Spectra 




10 3 10* 10 3 10* 



Multipole Moment (£) 

Fig. 2. — IR point source clustering power spectra from Planck (diamonds), BLAST (squares) and BLAST / ACT cross-correlations 
(triangles). The solid lines are the best-fit power law with scale dependence C^ lust oc £~ n - we find n = 1.25 ± 0.06. The frequency 
dependence of the clustering SED is described by a modified blackbody (see equation 3 and Table 1). Poisson shot- noise has been 
subtracted for each panel (shot- noise levels from our fitting are given in Table 3). For the combined fit the x 2 /d e g ree of freedom is 
132/122. The contribution to the total x 2 from the individual spectra is included in each panel. Color correction factors of 1.07 and 0.63 
have been applied to the BLAST 350 and 500 (im spectra so they are directly comparable with the Planck 857 and 545 GHz spectra, 
respectively, as described in the text. 
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TABLE 2 

IR CLUSTERING PREDICTIONS FROM OUR TEMPLATE 
t{t + l)C7f LUST /2vr (/1K 2 ) 



TABLE 3 

1-D MARGINALIZED SHOT-NOISE LEVELS 



Band 


v eB (GHz) 


ii 
z 


= 2000 


ii 

K. — 


= 3000 


Planck 100 GHz 


104 


0.49 


± 


0.07 


0.67 


± 


0.09 


Planck 143 GHz 


146 


2.9 


± 


0.3 


3.9 


± 


0.4 


Planck 217 GHz 


226 


15 


± 


4 


61 


± 


6 


ACT 148 GHz 


150 


3.2 


± 


0.4 


1.1 


± 


0.5 


ACT 218 GHz 


220 


36 


± 


3 


19 


± 


4 


SPT 95 GHz 


99 


0.38 


± 


0.05 


0.51 


± 


0.07 


SPT 150 GHz 


156 


4.1 


± 


0.4 


5.5 


± 


0.6 


SPT 220 GHz 


221 


37 


± 


3 


50 


± 


5 



ing component with index n = 1.2. Our predictions are 
in excellent agreement with these results. Figure 4 shows 
ACT and SPT data at 220 GHz with our clustering tem- 
plate predictions over-plotted, along with an IR shot- 
noise component and a ACDM CMB power spectrum. 

We also make predictions for the spectral index a 
of the clustered component. At the CMB frequencies 
the clustering SED can be approximated as a power 
law, I c {v) oc v a . however, the CMB bands do not lie 
strictly in the Rayleigh- Jeans limit of the modified black- 
body adopted in our model. Consequently the equivalent 
power-law slope at 150 GHz differs from that at 220 GHz: 
we find ai 50 = 3.78 ± 0.07 and a 22a = 3.56 ± 0.07. We 
also calculate an effective spectral index between 150 and 
220 GHz, ai5o-220 = 3.68±0.07. This result is consistent 
with ACT and SPT findings of 3.69±0.14 and 3.58±0.08, 
respectively. 

We have considered only the frequency dependence of 
the IR clustering component, not the shot-noise. Investi- 
gating shot-noise frequency dependence is more difficult, 
due to stronger flux cut dependence and the fact that the 
Planck data are limited to large scales; future Herschel, 
Planck, ACT, SPT, and cross-correlation data will help 
overcome these issues. 

4.1. Validating assumptions 

While our single- index power-law clustering model pro- 
vides a good fit to the data in terms of x 2 , in this section 
we attempt to further validate our assumptions regard- 
ing the scale and frequency dependence of the clustered 
power. 



Band 



Marginalized Value 
(Jy 2 sr- 1 ) 



Bethermin et al. (2011) 
(Jy 2 sr' 1 ) 



Flanck 857 GHz 


2200 


± 


2500 


5920 


± 


370 


Flanck 545 bnz 


1700 


± 


700 


1150 


± 


90 


Flanck 600 GHz 


210 


± 


60 


138 


± 


22 


r~) / „. „ „ 7„ ni 7 /~i IT 

Flanck 217 GHz 


16 


± 


6 


12 


± 


3 


DT A OT 1 OKO , 

BLAbl zoU f-tm 


11600 


± 


2300 


11600 


± 


2100 


RT A 3^0 nm 




4. 


1ZUU 


ouou 


4. 


luoU 


BLAST 500 v>m 


1300 


± 


390 


1680 


± 


480 


250 fj.m X 350 (im 


5700 


± 


1400 








250 fj.m X 500 £tm 


1830 


± 


690 








350 fim x 500 /^m 


820 


± 


190 








250 pirn x 218 GHz 


240 


± 


130 








250 pirn x 148 GHz 


90 


± 


70 








350 pirn x 218 GHz 


240 


± 


70 








350 fim x 148 GHz 


110 


± 


40 








500 (im x 218 GHz 


66 


± 


30 








500 ^tm X 148 GHz 


55 


± 


20 









4.1.2. Photometric calibration 

The marginalized values for the photometric calibra- 
tion parameters from our MCMC chains are shown in 
Table 4 along with the nominal uncertainties. We would 
expect the calibration parameters to be consistent with 
unity within the nominal uncertainty. If this were not 
the case it could indicate that the modified blackbody 
spectrum in equation (3) was not a suitable description 
for the frequency dependence of the clustering amplitude, 
however all are consistent with the nominal values. 



TABLE 4 

1-D MARGINALIZED CALIBRATION PARAMETERS 



Band 



Marginalized Value Nominal Uncertainty 



Planck 857 GHz 


1.05 


± 


0.06 


7% 


Planck 545 GHz 


0.99 


± 


0.04 


7% 


Planck 353 GHz 


1.01 


± 


0.02 


2% 


Planck 217 GHz 


1.00 


± 


0.02 


2% 


BLAST 250 um 


1.05 


± 


0.08 


9.5% 


BLAST 350 ura 


1.03 


± 


0.07 


8.7% 


BLAST 500 Atm 


1.02 


± 


0.07 


9.2% 


ACT 218 GHz 


1.03 


± 


0.07 


7% 


ACT 148 GHz 


1.00 


± 


0.02 


2% 



4.1.1. Shot-noise 

Table 3 shows the TD marginalized values of the shot- 
noise levels from our MCMC chains (Cp from equation 
1). It should be noted that, as expected, the Planck 
shot-noise levels and the clustering index n arc highly 
correlated. The shot-noise predictions for each auto- 
spectrum from the Bll model are given for comparison. 
Our values are consistent with the model within 1.5cr in 
all cases. We repeated the MCMC analysis described 
in Section 3.2 using the Bll model predictions as Gaus- 
sian priors on the Planck shot-noise levels and found a 
best-fit x 2 / d -°-f- = 139/122, an increase of Ax 2 = 7 
compared to the case with uniform priors. There was 
minimal change in the clustering template parameters 
(e.g., n — 1.24 ± 0.03 with the shot-noise prior). This 
suggests that uncertainty in the Planck shot-noise levels 
is not having a significant impact on our results. 



4.1.3. Planck beam uncertainty 

As stated in Section 2, we neglected any possible cor- 
relation in the Planck beam uncertainty across different 
angular scales. We re-ran the MCMC chain enforcing a 
100% correlation in the beam uncertainty in different £- 
bins for each Planck spectrum (introducing off-diagonal 
elements in the likelihood - equation 7) and found mini- 
mal change in the marginalized parameter values (< 0.3a 
in all cases). While the Planck beam uncertainties at 
I ~ 2000 are comparable to or larger than the statisti- 
cal uncertainties, the template parameters are primarily 
constrained by the large-scale Planck data (where the 
beam uncertainties are sub-dominant) and the BLAST 
data. 

4.1.4. Comparison with broken power law 
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Fig. 3. — Frequency dependence of IR point source clustering power at I = 3000. We assume this frequency dependence can be described 
by a modified blackbody: C^|qq (^) oc [u^B(u, T e g)] 2 . The solid lines show the best-fit and 1-ct uncertainties from our fitting; parameter 
values are given in Table 1. The dashed line is the SED of Gispcrt et al. (2000), which has been scaled for comparison. The rectangles 
show the Planck, BLAST, ACT and SPT bandpass filters F WHM values in the horizontal direction. The vertical extent of the rectangles 
shows the variation of the clustering power across each filter; in some cases this variation is a factor of 10 or more due to how steeply the 
clustering SED rises with frequency. The shaded rectangles show the bands that were used for the fitting in this paper. Spectra from the 
other bands are either not yet available or were not used, due to the risk of biasing results by assumptions regarding the separation of the 
CMB and other components from the IR contribution. 



We assumed that the scale dependence of the clustering 
power can be described using a single index n. In this 
section we consider fitting the scale dependence with a 
broken power law of the form 



a 



clust 



e~ ni if i < 4 

£-n 2 if £ > - 



b, 



(8) 



where m, n% and £b are free parameters. When we repeat 
the fit to the Planck, BLAST and BLAST x ACT spectra 
with this new scale dependence we find n\ = 1.24 ±0.06, 
n 2 = 1.2±0.2 and £ b = 1100±300 with a best X 2 /d.o.f. = 
131/120. This is an improvement of only A% 2 = 0.9 with 
two additional fitted parameters. For the fit with priors 
on the Planck shot- noise levels from the Bll model the 



improvement is only A\ 2 = 0.4. We conclude that the 
data do not show a significant preference for a broken 
p ower law. 

iKeisler et alj (|201lD adopted a template of the form 
given in equation (8) with n\ = 2.0, 712 = 1.2 and t\, = 
1500 to model the IR clustered power when extracting 
cosmological information from the small-scale SPT 150 
GHz CMB power spectrum. We find that this form is 
not a good fit to the Planck data, with x 2 /d.o.f. = 22/7 
when we fit to the Planck 217 GHz spectrum with n±, 
JI2 and £h fixed to the above values and using uniform 
priors on the clustering amplitude and shot-noise. 

4.1.5. Comparison with linear bias model 



Clustering template 



9 



HI 2 found that the BLAST and BLAST x ACT data 
arc well-fit by assuming the IR galaxy power spectrum 
in equation (5) is given by 

P g ai(M) = b 2 P DM (k,z) 7 (9) 

where b is the linear bias factor and Pdm is the linear 
dark matter power spectrum, and using the Bll predic- 
tions for the redshift-distribution of the flux, dS/dz. Pll 
also found that the Planck data from each band are well- 
fit by this model (also using the Bll dS/dz) if no prior is 
enforced on the shot-noise levels. The bias levels are not 
consistent; when fitting for a single, redshift-independent 
value of b, H12 found b = 5.0 ± 0.4 whereas Pll found 
b = 2.18 ± 0.11 for the Planck 545 GHz spectrum. 

The linear bias model is strongly rejected by the com- 
bined Planck and BLAST data even without any shot- 
noise prior. Fitting a linear bias model with single-value 
bias to the Planck 857 GHz and BLAST 350 /im spec- 
tra together, using the Bll dS/dz and multiplying the 
BLAST data points by a color correction factor of 1.07 
(see Section 4), yields x 2 /d.o.f. = 39/16. For a joint fit 
to the Planck 545 GHz and BLAST 500 /j,m spectra we 
find x 2 /d.o.f. = 43/16. This result is driven by the shape 
of the linear matter power spectrum rather than the 
choice of dS/dz: we repeated the fitting using the predic- 
tions of lMarsden et al.l(|2011l ) rather than Bll and found 
X 2 /d.o.f. of 36/16 and 42/16 for the 857 GHz / 350 fim 
and 545 GHz / 500 /im spectra respectively. Neither the 
Planck nor BLAST data alone were able to rule out the 
linear bias model without a shot- noise prior because o f 
the limited angular scales probed. IShirokoff et al.l (|2011[ ) 
similarly found that the linear bias model could not be 
ruled out using only SPT data from £ > 2000. 

4.2. Cross- correlations 

Table 5 shows the marginalized degree of clustering 
correlation, / CO rr, for each cross-spectrum included in 
our fitting, along with the separation Ais between the 
two bands. The BLAST x BLAST cross-spectra are 
consistent with 100% correlation. Our conclusions re- 
garding the degrees of correlation between the ACT and 
BLAST bands are limited by the data quality. A de- 
crease in correlation with increasing band separation 
would be consistent with the sources lying at a range 
of redshifts, with the higher-redshift sources being of 
great er relative importance at the longer wavelengths 
(e.g.. lHaiman fc Knoxj [2000l ) . however the current data 
are not of sufficient quality to confirm this. 

Two of the marginalized mean f co „ values lie more 
than lcr above unity, while none lie more than lcr be- 
low. Measuring / cor r > 1 may indicate that the angular 
scale dependence of the cross-spectra clustering power is 
not described by the same single-index power law as the 
auto-spectra clustering, since there is no physical expla- 
nation for a correlation in excess of 100%. This could be 
the case even with no worsening of the x 2 , due to the lim- 
ited angular scales probed by the BLAST x ACT data. 
To test that this is not having a significant effect on our 
clustering template, we repeated the MCMC fitting de- 
scribed in Section 3.2 using only the auto-spectrum data. 
We found that the values of n, /3 and Iq change by < 0. So- 
compared to the fit with the cross-spectra included. We 
conclude that, while the data may be hinting that the 



cross-spectrum clustering power has a different shape to 
the auto-spectrum clustering, this is not significantly bi- 
asing our template. 

Further submm-mm cross-correlation studies (e.g., 
Herschel x ACT / SPT) are clearly required to pro- 
vide more insight into the distribution of the sources 
with redshift, to constrain the angular scale dependence 
of the cross-spectrum clustering, and to investigate how 
the clustering shape changes with increasing band sepa- 
ration. 

TABLE 5 

Degree of clustering cross-correlation, / CO rr 



Band 1 Band 2 Au (GHz) Marginalized Value 



250 


(im 


350 fim 


340 


0.98 


± 


0.30 


250 


fim 


500 fim 


600 


1.23 


± 


0.31 


350 


fim 


500 fim 


260 


1.41 


± 


0.28 


250 


fim 


218 GHz 


982 


0.66 


± 


1.09 


250 


fim 


148 GHz 


1052 


0.30 


± 


1.83 


350 


fim 


218 GHz 


642 


0.64 


± 


0.56 


350 


fim 


148 GHz 


712 


0.39 


± 


1.05 


500 


fim 


218 GHz 


382 


1.82 


± 


0.47 


500 


fim 


148 GHz 


452 


0.79 


± 


0.87 



5. CONCLUSIONS 

We have found that a power-law model for the IR 
point source clustering is adequate to simultaneously 
fit Planck, BLAST and cross-correlated BLAST / ACT 
power spectrum data over a broad range of frequency 
(150 < v < 1200 GHz) and angular scale (multipole mo- 
ment 100 < £ < 9000). We find the clustering power 
varies with angular scale as l~ n with n = 1.25 ± 0.06 
and that the SED of the clustering can be described as a 
modified blackbody with cmissivity index j3 — 2.20±0.07 
and effective temperature T c ff = 9.7 K. Our work does 
not rely on any assumptions regarding the physical prop- 
erties of the IR sources (host halos, redshift distribution, 
etc.). 

As well as providing a simple template for use in CMB 
foreground subtraction, we have established that the 
Planck and BLAST / BLAST x ACT data sets appear 
compatible when bandpass filters, flux cut and calibra- 
tion are accounted for, as described in Sections 2 and 
3. We make predictions for the IR clustering power 
for the Planck, ACT and SPT CMB bands; our predic- 
tions for the ACT and SPT bands at around 150 and 
220 GHz are fully consistent with existing measurements 
(|Dunklev et alJl201Tl: [ShTrokoff et al.ll2^Ll . 

There is uncertainty in the Planck shot-noisc levels be- 
cause Planck does not probe small enough scales for the 
shot- noise to become dominant. We find reasonable con- 
sistency between the shot-noise levels from our fitting 
and the predict ions of the parametric IR galaxy evo- 
lution model of iBethermin et all (|2011[ ). We repeated 
our fitting using this model's predictions as priors on the 
Planck shot-noise levels and found there was minimal 
effect on the scale or frequency dependence of the clus- 
tering power. Number counts obtained via P{D) analysis 
(probability of deflection - modeling the probability dis- 
trib ution function of observed flux in each ma p pixel - 
e.g., iPatanchon eTall [200l iGlenn et al.l [2010h may be 
useful for constraining shot-noise levels in future work. 
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Fig. 4.— The angular power spectra at 220 GHz measured by ACT l IDas et al.H2011l ) and SPT l lShirokoff et al.|[201lT ), with a theoretical 
model for CMB and infrared point sources over-plotted. The clustered IR component is given by the power-law model fit to Planck, 
BLAST and BLAST X ACT data in this analysis (see equations f and 3 and Table 1 for parameter values). The upper and lower dashed 
lines correspond to l a error bounds. The lensed CMB spectrum is that of the ACDM model with parameters derived from WMAP 
(Komatsu ct al. 201f). An IR shot-noise component of size l)Cp/27r|^— 3000 = 78 fiK 2 (consistent with ACT and SPT measurements) 
has also been plotted. Additional radio source and kinetic SZ power arc subdominant at this frequency and are not included. 



Upcoming data from Herschel, Planck, ACT, SPT. 
and cross-correlations will allow us to look for any vari- 
ations in the clustering angular scale dependence with 
frequency. Understanding this scale dependence in the 
mm CMB bands is important not only for extracting 
unbiased estimates of cosmological parameters, but also 
for constraining other components of the measured spec- 
trum, such as the thermal and kinetic SZ effect, and any 
cross-component correlations, for instance SZ-IR. 

Over the range of angular scales considered we would 
expect to be probing both the linear and non-linear clus- 
tering regimes. This is supported by the fact that the 
combined Planck and BLAST data strongly reject a lin- 
ear bias clustering model (Section 4.1.5). The linear and 
non-linear components have different scale dependence, 
with the linear component dominating on large scales and 
the non-linear on small scales. The fact that our study 
has found that a power-law scale dependence provides 
a good fit to the Planck and BLAST IR point source 
clustering power is not inconsistent with this picture; in- 
stead it suggests that the sum of the linear and non- linear 
components is sufficiently close to a power law that the 
current data are unable to reveal any deviations. We 
can draw a comparison to measurements of the angular 



correlation function, w(9), of luminous red galaxies and 
galaxies in the main SDSS sample, which have revealed 
deviations from power law b ehavior (e.g., iZehavi et al.l 
12001 12001 iBlake et al.ll200l whereas earl ier studies of 
galax y clustering were unable to do so (e.g. . IZehavi et al.l 
I2002T ). We may expect similar deviations to be detected 
in future IR source power spectra, in particular with Her- 
schel /SPIRE data, which will yield higher quality data 
than BLAST out to smaller scales with its improved an- 
gular resolution. 

Th e correlation function of resolved sub-m m sources 
(e.g., iCoorav et all 120101 : iMaddox et all I2010T ) provides 
complementary information to the power spectra of un- 
resolved so urces considered in this work. Current data 
is limited; IGuo et aT] (|201lD measure £(r) cx r~ 7 with 
7^2, corresponding to Cf ust cx t~ x , for low-redshift 
(z < 0.5) Herschel- ATLAS galaxies, although their re- 
sult is consistent with ours within errors . The full Her- 
schel-ATLAS survey (|Eales et all I2010D will cover 30 
times more sky and provide far tighter constraints. Dis- 
crepancies between the Cf ust and £(r) or w{9) measure- 
ments would indicate differences in clustering properties 
of the bright sub-mm galaxies compared to the fainter 
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population, and thus both statistics are important for 
constraining models of galaxy clustering and evolution. 
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